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Summary: 

Three new iteration methods, namely the squared-operator method, the modified squared- 
operator method, and the power-conserving squared-operator method, for solitary waves 
in general scalar and vector nonlinear wave equations are proposed. These methods are 
based on iterating new differential equations whose linearization operators are squares of 
those for the original equations, together with acceleration techniques. The first two meth- 
ods keep the propagation constants fixed, while the third method keeps the powers (or 
other arbitrary functionals) of the solution fixed. It is proved that all these methods are 
guaranteed to converge to any solitary wave (either ground state or not) as long as the 
initial condition is sufficiently close to the corresponding exact solution, and the time step 
in the iteration schemes is below a certain threshold value. Furthermore, these schemes are 
fast-converging, highly accurate, and easy to implement. If the solitary wave exists only 
at isolated propagation constant values, the corresponding squared-operator methods are 
developed as well. These methods are applied to various solitary wave problems of physical 
interest, such as higher- gap vortex solitons in the two-dimensional nonlinear Schrodinger 
equations with periodic potentials, and isolated solitons in Ginzburg-Landau equations, 
and some new types of solitary wave solutions are obtained. It is also demonstrated that 
the modified squared-operator method delivers the best performance among the methods 
proposed in this article. 



1 Introduction 

Solitary waves play an important role in nonlinear wave equations. While such waves in 
some wave equations can be obtained analytically (such as in integrable equations), they 
defy analytical expressions in most other cases. Thus numerical computations of solitary 
waves is an important issue. In the past, a number of numerical methods have been de- 
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veloped for solitary waves. One of them is Newton's iteration method. This method can 
converge very fast. However, it often employs the finite-difference discretization, which has 
a low accuracy (compared to spectral or pseudospectral methods). In addition, in higher 
dimensions, time-efficient programming of this method becomes more cumbersome. Recent 
studies also showed that this method can suffer erratic failures due to small denominators 
[1]. Another method is the shooting method (see [2] for instance). This method works 
for all one-dimensional problems and higher-dimensional problems which can be reduced 
to the one-dimensional problem (by symmetry reduction, for instance). It is efficient and 
highly accurate. However, it fails completely for higher-dimensional problems which are not 
reducible to the one-dimensional problem. A third method is the nonlinear Rayleigh-Ritz 
iteration method (also called the self-consistency method) [3, 4], where one treats the non- 
linear eigenvalue problem as a linear one with a solution-dependent potential. The solitary 
wave is obtained by repeatedly solving the linear eigenvalue problem and normalizing the 
solution. This method can become cumbersome as well in high dimensions when the lin- 
ear eigenvalue problem becomes harder to solve. Two more methods are the Petviashvili 
method [5, 6, 7, 8] and the imaginary-time evolution method [9, 10, 11, 12, 13]. The con- 
vergence properties of the former method were studied in [6] for a class of equations with 
power nonlinearity, while those of the latter method were obtained in [13] for a much larger 
class of equations with arbitrary nonlinearity. Interestingly, it was shown that the conver- 
gence of the latter method is directly linked to the linear stability of the solitary wave if 
the solitary wave is sign-definite [13]. Both methods converge fast, are highly accurate, 
are insensitive to the number of dimensions, and their performances are comparable [13]. 
However, both methods diverge for many solitary waves (especially ones which cross zeros, 
i.e., excited states) [6, 13]. In recent years, some interesting yet challenging solitary wave 
problems arise in physical applications. Two notable examples are nonlinear light propaga- 
tion in multi-dimensional periodic and quasi-periodic media, and Bose-Einstein condensates 
loaded in multi-dimensional harmonic magnetic traps and periodic optical lattices. These 
problems arc not reducible to one-dimensional problems, so the shooting method can not 
be used. In addition, solitary waves in these problems often cross zeros (as is always the 
case in higher bandgaps), thus the Petviashvili method and the imaginary-time evolution 
method do not converge. Furthermore, the numerical stability analysis of such solutions 
require the solutions themselves to have high accuracy, which is often hard to achieve by the 
Newton's method or the self-consistency method. Thus new numerical schemes are called 
upon. These schemes should be time-efficient, highly accurate, insensitive to the number 
of dimensions, and capable of computing all types of solitary waves in any scalar or vector 
nonlinear wave equations. Equally importantly, these schemes should be easy to implement. 
None of the previous methods meets all these requirements. 
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In this paper, we propose several classes of iteration methods for solitary waves which do 
meet all the requirements described above. These methods are based on two key ideas. 
The first key idea is to iterate a modified equation whose linearization operator is a square 
of that for the original equation. This idea (in a different form) was first presented in 
[9]. Specifically, those authors proposed to obtain excited-state solitary waves as stationary 
points of a certain functional that is different from the usual Lagrangian. In the present 
study, we will show that this procedure is equivalent to using the aforementioned squared 
operator in the iteration equation. We further show that this operator structure guarantees 
that the proposed methods converge for all types of solitary waves in any nonlinear wave 
equations. These iteration methods are compatible with the pseudo-spectral method for the 
computation of spatial derivatives, thus they are highly accurate (with errors that decrease 
exponentially with the decrease of the spatial grid size), and can handle multi-dimensional 
problems with little change in the programming. The second key idea is to introduce an 
acceleration operator to the scheme [9, 13]. The acceleration operator speeds up the con- 
vergence of these iteration methods drastically, hence making them highly time-efficient. 
Based on these two key ideas, we propose two powerful new iteration methods which we 
call the squared-operator method and the power (or arbitrary quantity)-conserving squared- 
operator method. The former method specifies the propagation constant of the solitary 
wave, while the latter method specifies the power or other arbitrary functional of the solu- 
tion. Both methods are shown to converge to any solitary wave if the time-step parameter in 
the methods is below a certain threshold value, and the initial condition is sufficiently close 
to the exact solution. Beyond these two ideas, we also employ an eigenmode-elimination 
technique [14] which speeds up the convergence of iterations even further. The resulting 
numerical method will be called the modified squared-operator method in the text. All 
these schemes pertain to solitary waves which exist at continuous propagation constants. 
In certain nonlinear wave problems (especially of dissipative nature), however, solitary 
waves exist at isolated propagation constants. By extending the above ideas, we construct 
the squared-operator method and the modified squared-operator method for isolated soli- 
tary waves as well. By applying these new schemes to various challenging solitary wave 
problems, we demonstrate that the modified squared-operator method gives the best per- 
formance, followed by the squared-operator method, then followed by the power-conserving 
squared-operator method. In the end, we construct a whole family of iteration schemes 
which share similar convergence properties and contain the presented schemes as special 
cases. 

This paper is structured as follows. In Section 2, we present the squared-operator method 
and show that it converges for all types of solitary waves in arbitrary vector equations. We 
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also demonstrate it on the nonlinear Schrodinger equation where we carry out the explicit 
convergence analysis and determine optimal parameters in the acceleration operator. In 
Section 3, we present the modified squared-operator method and show that this method not 
only converges for all types of solitary waves in arbitrary equations, but also converge faster 
than the squared-operator method. In Section 4, we present power-conserving versions of the 
squared-operator method and prove their universal convergence properties. In Section 5, we 
derive the squared-operator method and the modified squared-operator method for isolated 
solitary waves. In Section 6, we demonstrate all these methods on various challenging 
examples, including vortex-array solitons in higher bandgaps and solitary waves in coherent 
three- wave interactions, and show that the modified squared-operator method delivers the 
best performance. In the Appendix, we present whole families of squared-operator-like 
methods, and show that the methods described in the main text are the leading members 
in these families in terms of convergence speeds. 



2 The squared-operator method for solitary waves in general 
nonlinear wave equations 

Consider solitary waves in a general real- valued coupled nonlinear wave system in arbitrary 
spatial dimensions, which can be written in the following form 

Lou(x) = 0. (1) 

Here x is a vector spatial variable, u(x) is a real-valued vector solitary wave solution ad- 
mitted by Eq. (1), and u ^ as |x| oc. Note that for complex-valued solitary waves, 
the equation can be rewritten in the above form with u containing the real and imaginary 
parts of the complex solution. Let Li denote the linearized operator of Eq. (1) around the 
solution u, i.e., 

Lo (u + u) = Liu + 0(u2), u < 1, (2) 

(■)^ denote the Hermitian of the underlying quantity, and M be a real- valued positive- 
definite Hermitian operator which we call the acceleration operator. The idea of our method 
is to numerically integrate the time-dependent equation 

ut = -M-^l|m-^Lou (3) 

rather than Uj = ibM~^LoU. The reason is to guarantee the convergence of numerical 
integrations, as we will demonstrate below. The operator M is introduced to speed up the 
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convergence, in the same spirit as the pre-conditioning technique in solving systems of hnear 
algebraic equations. Using the simplest time-stepping method for Eq. (3) — the forward 
Euler method, the iteration method we propose for computing solitary waves u is 



M 



At. 



u=u„ 



(4) 



It will be shown below that this method is universally convergent as long as the time step 
At is below a certain threshold, and this universal convergence stems from the fact that 
the iteration operator for the error function is — M~'''L|M~^Li, or "square" of the operator 
M^-'^Li. Thus we call scheme (4) the squared-operator method (SOM). If M is taken as 
the identity operator (no acceleration), then SOM (4) reduces to 



Un+l = Un — 



lJLqU 



u=u„ 



At, 



(5) 



which has a simpler appearance but converges very slowly. 



It should be noted that even though more complicated time-stepping methods (such as 
Runge-Kutta methods) can also be used to integrate Eq. (3) , they are actually less efficient 
than the forward Euler method (4), because the extra computations in them outweighs the 
benefits they may have. Implicit methods are even less attractive. The reason is that due 
to the acceleration operator M, the time steps At allowed by explicit methods such as (4) 
for numerical stability (or convergence) are not small, thus the need for implicit methods 
vanishes. 

Scheme (4) is actually one of the many SOM's one can construct using the same squared- 
operator idea. Indeed, in the Appendix, we will present a whole family of SOM's, of which 
(4) is a particular member. We will show there that scheme (4) is the leading member in 
that family in terms of convergence speeds. 

Let us now remark on the relation of these methods to the functional minimization method 
for Hamiltonian equations proposed in [9] , which has the following form: 

6 



SuJ 



LquII dx. 



(6) 



Here 6/Su represents the functional (Frechet) derivative, and ||...|| denotes the L2-norm. 
Upon taking this functional derivative and noticing that for Hamiltonian equations, lJ = 



Li, one recovers Eq. (3) with M 
the functional equation 



1. The accelerated equation (3) follows similarly from 



ut = y ||m-i/2LoM- 



1/2 



u 



(7) 
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where u = M^/^u. 

To formulate the convergence theorem for the SOM (4), we introduce the following assump- 
tion. 

Assumption 1 Let w(x; ci, C2, . . . , Cg) he the general solution of Eq. (1), where (ci, C2, . . . , c^) 
are free real parameters. Then we assume that the only eigenf unctions in the kernel of Li 
are the s invariance modes Ucj,l < j < s, where Uc^ = du/dcj. 

Remark 2.1 One example of these invariance modes is the translational-invariance mode. 
Suppose Eq. (1) is translationally invariant along the direction xi, i.e., u{xi+ci,X2, ■■■ , xn) 
is a solution for any value of ci. Then the translational-invariance mode = Ua;^ is in 
the kernel of Li. Another example of the invariance modes pertains to arbitrary phases of 
complex- valued solutions and can be found, e.g., in Example 6.1 below. 

Remark 2.2 Assumption 1 holds in the generic case. However, in certain non-generic cases, 
it does break down (see Fig. 7 of [15] for an example). 

Under Assumption 1, we have the following theorem. 



Theorem 1 Let Assumption 1 he valid, and define 

_ 2 

^tmax = — T ) (8) 

^'-min 

where Kmin is the minimum eigenvalue of operator 

C = -M-^lIm-^Li. (9) 

Then SOM (4) is guaranteed to converge to the solitary wave m(x) of Eq. (1) if At < At max 
and the initial condition is close to w(x). If At > Atmax, then SOM (4) diverges from w(x). 



Proof. We use the linearization technique to analyze SOM (4) and prove the theorem. Let 

u„ = u -I- u„, u„(x) < 1, (10) 

where u„(x) is the error. When Eq. (10) is substituted into SOM (4) and only terms of 
0(u n) are retained, we find that the error satisfies the following linear equation: 

Un+l = {1 + At jC) Un, (11) 
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where L is defined in Eq. (9). Note that since £ = M-i/2£/,Mi/^ where 

C,n = - (M-V^LiM-V^y (m-V2LiM-V2) 

is Hermitian and scmi-negativc-dcfinite, thus all eigenvalues of L are real and non-positive. 
In addition, all cigenfunctions of £ form a complete set in the square-integrable functional 
space, hence can be expanded over £'s eigenf unctions. Consequently, if At > At^j^, the 
eigenmode with eigenvalue A^j„ in the error will grow exponentially due to l+A^j^At < — 1, 
i.e., SOM (4) diverges. On the other hand, if At < Atmim then no eigenmode in the error 
can grow. In fact, all eigenmodes decay with iterations except those in the kernel of C But 
according to Assumption 1, eigenfunctions in the kernel of C are all invariance modes Uc^ 
which only lead to another solution with slightly shifted constants cj and do not affect the 
convergence of iterations. Thus Theorem 1 is proved. 

Remark 2.3 The role of the acceleration operator M is to make A^j„ of £. bounded 
(without M, Amin = — oo in general); see, e.g.. Example 2.1 below and [13]. As shown after 
Remark 2.4 below, this leads to faster convergence of SOM (4). 

Remark 2.4 In computer implementations of SOM (4), spatial discretization is used. In 
such a case, if we require the discretization M*-^-* of the positive-definite Hermitian operator 
M to remain positive-definite and Hermitian (which is needed to show the non-positiveness 
of eigenvalues of the discretized operator JC^^^), then Theorem 1 and its proof can be 
readily extended to the spatially-discretized SOM scheme (4), except that Amin in Eq- (8) 
becomes the smallest eigenvalue of the discrete operator C^^^ now. Unlike discretizations of 
some other schemes such as the Petviashvili method [6] and the accelerated imaginary-time 
evolution method [13], discretizations of SOM (4) always converge regardless of whether 
the discretized solitary wave is site-centered (on-site) or inter-site-centered (off-site). Note 
that under discretization, the kernel of £^^^ may become smaller than that of £ due to the 
breakdown of translational invariances, but this does not affect the extension of Theorem 
1 and its proof to the discretized SOM scheme at all. 

In the application of SOM (4), a practical question is how to choose At within the upper 
bound of Atjnax so that convergence is the fastest. To answer this question, let us define 
the convergence factor R of SOM (4) as 

ii: = max{|l + AAt|}, (12) 

where A is any non-zero eigenvalue of £.. Then the error u„ of the iterated solution decays 
as i?", where n is the number of iterations. Smaller R gives faster convergence. Clearly, 

R = max{|l + A^i„At|, |l + A^„^At|}, (13) 
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where A^^j. is the largest non-zero eigenvalue of operator L. From this equation, we see 
that the smallest R (fastest convergence) occurs at the time step 

for which the corresponding convergence factor is 

'■^mm ~ '■^max 

Another practical question which arises in the implementation of SOM (4) is the choice 
of the acceleration operator M. This operator should be chosen such that M is easily 
invertible. In addition, it should be chosen such that the convergence is maximized or 
nearly maximized. A simple but often effective choice is to take M in the form of the linear 
part of the original equation (1), to which one adds a constant to make it positive definite. 

To demonstrate the effect of M on convergence speeds, below we consider the familiar 
nonlinear Schrodinger (NLS) equation for which the convergence analysis of SOM can be 
done explicitly. 

Example 2.1 Consider the NLS equation in one spatial dimension: 

Uxx + = nu. (16) 

Without loss of generality, we take fi = 1. For this fi value, Eq. (16) has a soliton solution 

u{x) = V^sechx. (17) 

We take the acceleration operator M to be in the form of the linear part of Eq. (16), i.e., 

M = c-d^^, c>0, (18) 

which is easily invertible using the Fourier transform. Then the eigenvalue equation for 
operator M~^Li, 

M-^Li^ = XtP, (19) 

can be rewritten in the explicit form: 

^xx - Yipy + Yqr^sech^xV' = o. (20) 



8 




The continuous spectrum of this equation can be obtained by requiring the eigenfunction 
to be oscillatory at |a;| = oo, and this continuous spectrum is found to be 

AG for c> 1; 

(21) 

AG [-i,-l), for c< 1. 
The discrete eigenvalues of (20) satisfy the equation [16, 17] 

/l + Yq^-(2i + l) , j = 0,l,2,..., (22) 

where the right hand side should be non-negative. These discrete eigenvalues are plotted in 
Fig. 1(a). Note that A = is always a discrete eigenvalue with j = 1. This zero eigenvalue 
is due to the translational invariance of the NLS soliton (which leads to LiUx = 0). The 
eigenvalues Aq and A2 (for j = and 2) can be found to have the following expressions: 

Ao(c) = 2(0 + 5f _ 

2c2 + 9c + 1 + V25c2 + 118c+l 

A2(c) = ^'^'^P' - 1. (24) 

2c2 - 3c + 85 + 5V25c2 - 26c + 145 

From the above spectrum of operator M~^Li, we can easily obtain the spectrum of the 
iteration operator jC = —M~^l\m~^Li, whose eigenvalues A are related to the eigenvalues 
A of M~^Li as A = — A^, since in this case l\ = Li. Hence 

^maxic) = -Ai(c), A^j„(c) = min |-Ao(c), -l| . (25) 

These eigenvalues are plotted in Fig. 1(b). Based on these eigenvalues, we can calculate 
the convergence factor i?*(c) from formula (15). This R*{c) is shown in Fig. 1(c). It is 
seen that when Ao(c) = 1, i.e., c = Copt = 6 — \/T3 ~ 2.4, i?*(c) reaches its minimum 
value Rmin ~ 0.80. At c = Cgpt, R*{c) is not differentiable, because Amm(c) equals — Ao(c) 
for c < Copt aiid ~1 (the edge of the continuous spectrum) for c > Copt- For c = Cgpt, 
the dependence of the convergence factor R on At is shown in Fig. 1(d). At At^ 1.80 
from formula (14), R reaches its minimum value Rmin given above. Without acceleration 
(M = 1), this value of R would be very close to 1 with discretizations (since A^m is large 
negative, see [13]), and be exactly equal to 1 without discretizations (since Amin = —00). 
Prom the above explicit analysis, we see that the choice of M affects the convergence speed 
of SOM (4) significantly. 
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3 The modified squared-operator method for sohtary waves 
in general nonhnear wave equations 



The above SOM (4), even with a sensible choice of the acceleration operator M, can still be 
quite slow for certain problems (see the example in Fig. 4). In this section, we employ an 
additional technique which can speed up the convergence of SOM iterations even further. 
This technique is called eigenmode elimination, and was originally proposed in [14] for a 
non-squared-operator scheme. When this eigenmode-elimination technique is incorporated 
into SOM, the resulting method, which we call as the modified squared-operator method 
(MSOM), will converge faster than SOM (4). 

The MSOM we propose is 

M-iLlM-iLoU-an(G„, l1m-1Lou)gJ _ At, (1) 



Un+l — U„ 

where 



u=u„ 



_ 1 1 

" (MG„,G„) " (LiG„, M-iLiG„)At' ^ ^ 

Gn is a function the user can specify, and the inner product is the standard one in the 
square-integrable functional space: 

(Fi,F2) = r f1 -Fsdx. (3) 

J —oo 

If (Fi,F2) = 0, Fi and F2 are said to be orthogonal to each other. Two simple choices for 
Gn can be 

Gn = u„, (4) 

and 

G„ = e„ = u„ - u„_i. (5) 

The motivation for MSOM (1) can be explained briefly as follows [14]. Consider SOM (4), 
and denote the slowest-decaying eigenmode of C in the error as G(x) with eigenvalue A^. 
Note that according to Eqs. (12) and (13), Ag = Amin or Amax- Our idea is to construct a 
modified linearized iteration operator for the error so that eigenmode G(x) decays quickly, 
while the decay rates of other eigenmodes of JC remain the same. If so, then this modified 
iteration scheme would converge faster than the original SOM. For this purpose, consider 
the modified linearized iteration operator 

Cm^ = jC^- a{MG, jC^)G, (6) 
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where a is a constant. Since G(x) is an eigenfunction oi C, and recalling that eigenfunctions 
of this C are orthogonal to each other under the M- weighted inner product, we readily see 
that this modified iteration operator and the original one have identical eigenfunctions. 
Their eigenvalues arc identical too except the one for eigenmode G(x). The eigenvalue of 
this eigenmode changes from of £ to (1 — q;(MG, G))A3 of Cm- Then if we choose a so 
that the decay factor of this eigenmode is zero, i.e., 

l + (l-a(MG,G))A,At = 0, (7) 

or equivalently, 

then eigenmode G(x), which decays the slowest in SOM, becomes decaying the fastest (in 
fact, this mode is eliminated from the error after a single iteration), while decay rates of the 
other eigenmodes in the error are unaffected. Thus convergence is improved. In practical 
situations, the slowest-decaying mode G(x) and its eigenvalue A^ in SOM are not known. 
In such cases, if we can choose G„ which resembles the slowest-decaying eigenmode of 
SOM, then the corresponding eigenvalue A^ can be approximated as the Rayleigh quotient 
Kg ~ (MG„, >CGn)/(MGn, Gn). Substituting this approximation into (8), we get a as 
given by formula (2). Corresponding to the modified linearization operator (6), the modified 
iteration method is then MSOM (1). 

The above derivation assumed that function G„ is equal or close to the slowest decaying 
eigenmode of SOM (4). In practical implementations of this method, one does not know 
priori if the selected G„ meets this criterion. This then may put the effectiveness of this 
method in question. One may also ask if this modified method can converge at all even 
with small time steps. Fortunately, the convergence of MSOM is insensitive to the choice of 
function G„, at least when certain mild conditions are satisfied. Specifically, let us expand 
function G„(u„) as 

G„ = Go + O(ii„), u„<l, (9) 

where Go is G„'s leading-order term, and u„ is defined by Eq. (10), then we have the 
following theorem. 



Theorem 2 Let Assumption 1 he valid, and Li Gq ^ 0, then if 



( 



At < AtM = min 



{MGM) "^'"7 



(10) 
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where C is given in Eq. ( 9), M is Hermitian and positive definite, and A^^^ is C 's minimum 
eigenvalue, then MSOM (1) is guaranteed to converge to the solitary wave u(x) of Eq. (1) 
if the initial condition is close to tx(x) . 

Proof. We again use the linearization technique to analyze MSOM and prove Theorem 
2. Substituting (10) into MSOM (1), we find that the error satisfies the linear iteration 
equation 

Un+i = (1 + At Cm) Un, (11) 

where 

^M* = C^- ao(MGo, i:*)Go, (12) 

and 

_ 1 1 
"° " (MGo,Go) ^ (MGo, >CGo)Ar ^^^^ 
Note that since M is Hermitian and positive-definite, ao is real. 

Below, we show that all eigenvalues of Cm are real and non-positive, the kernel of Cm 
is the same as that of Li, Cm has no square- integrable generalized eigenfunctions at zero 
eigenvalue, and under the time-step restriction (10), |1 + A^Atl < 1 for all non-zero eigen- 
values Am of Cm- Then under the assumptions of Theorem 2, MSOM (1) will converge, 
and Theorem 2 will then be proved. 

To proceed, we first write out the eigenvalue problem for operator Cm, which is 

- ao(MGo, C-^)Go = Am^-. (14) 

Taking the inner product between this equation and M£^', and recalling the form (9) of 
C and that ao is real and M Hermitian, we see that eigenvalues Am are all real. Next, we 
will analyze the eigenvalue problem (14) by expanding into eigenfunctions of operator 
>C, a technique which has been used before on other eigenvalue problems [18, 19]. First, 
notice from the proof of Theorem 1 that eigenvalues of C are all real and non-positive, and 
eigenfunctions of C form a complete set. Let the discrete and continuous eigenmodes of C 
be 

^VitW = Ajk'0fe(x), Afe < 0, A; = 1, . . . , m, (15) 

>CV(x; A) = AV'(x; A), A G / < 0, (16) 

where m is the number of Cs discrete eigenvalues, / is Cs continuous spectrum, and the 
orthogonality conditions among these eigenfunctions are 

{M^i,^j)=5i,j, (17) 
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and 

(MV^(x;A),V^(x;A'))=<5(A-A'). 
Then we expand Go and ^(x; Am) as 



Go = X! CkM^) + / c(A)V'(x; A)dA, 
fc=i 

m „ 

^'(x;Am) = V6fc^fc(x) + / 6(A)?/;(x; A)dA. 



(18) 



(19) 



(20) 



Substituting Eqs. (19) - (20) into the eigenvalue problem (14) and using the orthogonality 
conditions (17) - (18), one obtains 



aoCfc(MGo, 



6(A) 



aoc(A)(MGo, C^) 



Notice that 



Afc-AM A -Am 

(MGo, = Y.Akbk4+ / A6(A)c(A)*(iA 

k=l •'^ 

ao(MGo, m if,-^^+ f 



(21) 



A|c(A)|^ 



A- A 



dA 



thus we get 
where 



(MGo, -C^r) • Q(Am) = 0, 



A 1^ |2 

q{Am) = E -^^^ + 



Afc - Am Ji A- Am 
Then the discrete eigenmodes of jCm are such that either 

(MGo, JC^) = 0, 



(22) 

(23) 
(24) 

(25) 
(26) 

The continuous eigenvalues of Cm are the same as those of C since Cm ^ C as |x| oo. 

We first consider eigenvalues of Cm where condition (25) holds. In this case, Eq. (14) 
becomes the eigenvalue equation for C, thus these eigenvalues of Cm are also the eigenvalues 
of C. As a result, 

Amin < Am < 0. (27) 



or 



Q(A 



M) 



0. 
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Next, we consider eigenvalues of Cm which satisfy condition (26). For these eigenvalues, we 
will show that the following inequality holds: 

min (A^j„, 

ao + A^in) < Am < 0, (28) 

where 

ao = -ao(MGo,/:Go) = - (mGo,Go) " At" ^^'^ 

To prove the right half of this inequality, recall that Aj. < and A G / < 0, hence for any 
real number Km > 0, 



^ Ajfc - Am JiA-Am 7/ 



(30) 



On the other hand, 

J_ ^ (MGq, Go) ^^^^ 

an 1 -L (MGq, Go) 
^ + (MGo,£Go>At 

Since LiGq 7^ by assumption and M is Hermitian and positive definite, then (MGq, JCGq) < 
0. In addition, (MGq, Go) > 0. Thus 

— > (MGo,Go) or — < 0. (32) 
ao "0 

In view of Eqs. (24), (30) and (32), Q{Am) can only have negative roots. Thus, the right 
half of inequality (28) holds. 

In the following, we prove the left half of inequality (28). Notice that 

. (MV>, £^) 
^-^" = T (MV^,^) ' 

thus ao + Kjnin < in view of Eq. (29) and At > 0. Let us rewrite 1/ao as 
1 (MGo, CGo) 1 



ao ao ao 

When this expression is substituted into Eq. (24), we get 
1 



E^fe|cifcl' + y/|c(A)|2dAj . (34) 



g(AM) 

ao 



(Afc + ao- AM)Afc|cfcp /• (A + ao- Am)A|c(A)P ^^^ 



Afe - Am ji a - Am 



(35) 
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When Am < min (Aj„j„, ctQ + -^mm)) all terms inside the square bracket of Eq. (35) are 
negative, thus Q{Am) can not be zero. As a result, the left half of inequality (28) holds. 

Due to the two inequalities (27) and (28) on eigenvalues Am for the two cases (25) and 
(26), we see that all eigenvalues of Cm are non-positive, and the kernel of jCm is the same 
as the kernel of Li. In addition, the convergence condition |1 + A^ At| < 1 for non-zero 
eigenvalues Am of Cm will be satisfied if 



At < , (36) 



and 



At<-- . (37) 

Since (Xq + A^j„ < and due to Eq. (29), the second inequality (37) is equivalent to 

< (MG,.Go) ^ • (^«) 
(MGo,Co> ~ '"^^ 

Together with Eq. (36), we find that when At satisfies the restriction (10), |1 -|- AMAt| < 1 
for all non-zero eigenvalues Am of Cm- 

Lastly, using similar techniques as employed above, we can readily show that Cm has no 
generalized eigenfunctions at zero eigenvalue, i.e., equation 

CmF = Uc^, l<j<s (39) 

has no square-integrable solutions F. Here Uc^. is in the kernel of Cm, i-e. the kernel of Li 

(see above and Assumption 1). This rules out the possibility of linear (secular) growth of 
zero-eigenmodes Ucj in the iteration of the error function u„. This concludes the proof of 
Theorem 2. 

To illustrate the faster convergence of MSOM, we apply MSOM with the choice of (4) to 
the NLS equation (16), where explicit convergence analysis can be carried out. 

Example 3.1 Consider MSOM (1), (4) applied to the solitary wave (17) in the NLS 
equation (16). For simphcity, we take M as in (18) with c = /x = 1. In this case, by 
inserting c = 1 into Eq. (22), we find that the eigenvalues of operator M'^Li are 



24 

(2j +"3p - 1 



= TW-T^;^^^ - J = 0,1,2,... (40) 
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Hence Aq = 2, Ai = 0, A2 = —0.5, . . . , Aqo = —1- Eigenvalues of operator L are 

A = -A2, j = 0,1,2,... (41) 

Notice that M~^L\u = 2u, hence Go = u is an eigenfunction of JC. According to the 
discussions at the beginning of this section, we find that eigenvalues of jCm are 

Am,o = (42) 

AM,j = -Xl J = 1,2,... (43) 
Hence the convergence factor Rm for MSOM (1) is 

Rni^t) = max{|l - XlAt\,\l - A^At|} = max{|l - 0.25At|, |1 - At|} . (44) 

It is noted that the zero eigenvalue Am,i corresponds to a translational-invariance mode and 

docs not affect the convergence analysis. From Eq. (44), we see that MSOM (1) converges 
if At < Atmax = 2. The fastest convergence is achieved at 

and the corresponding convergence factor is 

This convergence factor is substantially lower than that of SOM (4) (see Fig. 1), thus 

MSOM (1), (4) converges much faster. We note in passing that this convergence factor (46) 
for MSOM is close to those of the optimal Petviashvili method and the optimally accelerated 
imaginary-time evolution method, which can be calculated to be 0.5 and 0.51 respectively 
[13]. 

Remark 3.1 Unlike Atmax in Theorem 1 for SOM, AtM in the time-step restriction (10) 
is not a sharp bound for convergence of MSOM (1). In practice. At can often be taken 
larger than AtM, and MSOM still converges. This can be seen clearly in Example 3.1, 
where A^m = —4, Go = u, {MGq, CGq) / {MGo,Go) = 4, hence AtM = |- But according 
to the above explicit calculations, the sharp bound on the time step is At^ax = 2, much 
larger than AtM- Thus, restriction (10) is sufficient, but not necessary, for the convergence 
of MSOM (1). 
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Of the two simple choices for G„ given in (4) and (5), Gq = u for the first choice (4), 
and Go = for the second choice (5). Thus Theorem 2 apphes to choice (4), but not to 
choice (5) since LiGq = there. The convergence analysis for MSOM with choice (5) is 
more complicated since the corresponding linearized iteration equation for the error is no 
longer a one-step iteration, but rather a two-step iteration. However, choice (5) has certain 
intuitive appeal, as the difference function — u„_i often resembles the slowest decaying 
mode in the error [14], thus it makes sense to choose (5). Our numerical experiments have 
shown that indeed, MSOM with (5) not only converges (when the time step is below a 
certain threshold), but also converges much faster than the choice (4) in almost all cases. 
In fact, we will see from various examples in Sec. 6 that MSOM (1) with choice (5) often 
gives the fastest convergence among schemes proposed in this paper, especially for solitary 
waves with complicated profiles. 



4 Power-conserving and arbitrary-quantity-conserving squared- 
operator iteration methods 



In some applications, solitary waves are sought with a pre-specified power (or other quantity 
such as energy) rather than the propagation constant. For instance, in Bosc-Einstcin con- 
densation, the number of atoms is often specified, and a solitary wave with that number of 
atoms (i.e. power in our notations) needs to be computed. In principle, we can still use the 
aforementioned SOM or MSOM to compute the wave by first continuously varying propaga- 
tion constants to obtain the power curve, then determining the propagation constant with 
the pre-specified power, then finally computing the wave again with SOM/MSOM. But that 
is clearly awkward and not time efficient. A much better way is to design numerical schemes 
which automatically converge to the solution with the pre-specified quantity. Another ex- 
ample is the linear eigenvalue problem. In this case, the eigenvalues are unknown, thus 
SOM/MSOM clearly does not apply. If one views the eigenvalues as propagation constants, 
and requires the eigenfunctions to have a fixed power (norm), then this linear eigenvalue 
problem becomes the same as a solitary-wave problem with a pre-specified power. To treat 
this type of problems, one can use the imaginary-time evolution method [9, 10, 11, 12, 13], 
where the solution is normalized at every step to keep the pre-specified power. However, 
the problem with the imaginary-time evolution method is that it often diverges when the 
solution crosses zero [13]. Thus new numerical schemes which conserve the power but also 
have guaranteed convergence need to be constructed. 
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In this section, we propose power- (or any other quantity-) conserving squared operator 
methods, where one specifies the power (or any other quantity) instead of the propagation 
constant. These methods are in the same spirit as the one presented in [9], but we go much 
further here. First, we will rigorously prove that these new power-conserving methods 
converge to any solitary wave as long as the time step is below a certain threshold, and 
the initial condition is close to that solution — just like SOM/MSOM. This has never been 
done before. We should point out that this guaranteed convergence for all solitary waves 
is nontrivial; it is certainly not automatic for an iteration method even if the "squared 
operator" idea has been incorporated. This guaranteed convergence is possible only when 
the updating formula for the propagation constants is compatible with the solution updating 
scheme, in addition to the " squared operator" idea. Second, the methods we will propose use 
a different acceleration than the acceleration of [9] , and hence, as we show in the Appendix, 
have faster convergence. Thirdly, our methods apply to all types of equations and arbitrary 
forms of conserved quantities, more general than the method of Ref. [9]. 

For the ease of presentation, we first consider three special (yet large) classes of equations 
and construct their power-conserving squared-operator methods (PCSOM). Then we con- 
sider the most general case where the wave equations and the pre-specified quantities of 
the solution are both arbitrary, and construct their quantity-conserving squared-operator 
methods (QCSOM). 

4.1 The power-conserving squared-operator method for equations with a 
single propagation constant 

Consider equation (1) of the form 



Here u(x; fi) is a real vector solitary wave, and /x is a real scalar constant (usually called 
the propagation constant). Solution u is assumed to exist for continuous ranges of /x values. 
Equations of the form (1) include all scalar equations as well as certain vector equations 
such as vortex equations (2)- (3) in Example 6.1 and the second-harmonic generation system 
(9)-(10) in Example 6.2. Define the power of solution u(x;/x) as 



then we are interested in finding a solution u(x; /x) whose power has a pre-specified value 
P. Combining ideas of SOM (4) and the imaginary-time evolution method [9, 10, 13], we 



Lqu = Loqu — /iu = 0. 



(1) 



^(m) = (u(x;/x),u(x;/x)) 



(2) 
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propose the following new power-conserving squared-operator method (PCSOM): 



Un+l 



Un+l, 



where 



u„+i = u„ - M ^ l\m ^Lqu - 7U 

L Ju=u„,/U=/i„ 

M is a positive-definite and Hermitian acceleration operator, and 
(u, M-^lIm-ILou) (u, M-^Loou) 



At, 



(u, M-iu) 



(u, M-iu) 



(3) 
(4) 

(5) 



The convergence property of this PCSOM is similar to that of the SOM, and is summarized 
in the following theorem. 

Theorem 3 Let Assumption 1 he valid, u(x) he orthogonal to the kernel of hi, P'{fJ,) = 

dP{^)/d^ 7^ 0, and define Atmax by Eq. (8), where A„iin now is the minimum eigenvalue 
of the iteration operator Cpc defined in Eq. (7) below. Then when At < At^axj PCSOM 
(3)- (5) is guaranteed to converge to the solitary wave u(x) of Eq. (1) with power P if the 
initial condition is close to u(x). When At > Atmax, PCSOM (3)- (5) diverges. 

Proof. As before, we use the linearization method to analyze PCSOM (3)-(5). Substituting 
Eq. (10) into the scheme and linearizing, we find that the iteration equation for the error 

Un is 

Un+l = Un + C,PC Un At, (6) 



where operator Cpc is 



and 



Cpc'^ = M-^ (l* - 7u) , 
= -L|m-i (Li* - I3u) , 



7 



(u, M-^L*) 



(u, M-^Li*) 



(u, M-iu) ' (u, M-iu) 

In addition, due to the power normalization (3), we also have the constraint: 



(7) 
(8) 

(9) 



(Un, u) = 0. 



(10) 
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Thus it is sufficient to consider eigenfunctions ^'(x) of operator Cpc in the restricted space 
S: 

5={*(x): (*,u)=0}. (11) 



First, we rewrite the operator Cpc as Cpc = M~^/^>Cpc//M^/-^, which defined operator 
J^PCH- Then Cpc and CpcH are similar, hence having the same eigenvalues. Since we only 
need to consider eigenfunctions of Cpc which are orthogonal to u, it follows that we only 
need to consider eigenfunctions of CpcH which are orthogonal to M^^/^u. In this space, it 
is easy to check that operator CpcH is Hermitian. In addition, using the Cauchy-Schwartz 
inequality, we can verify that CpcH is also semi-negative definite. Thus, all eigenvalues 
of Cpc are real and non-positive, and the eigenfunctions of Cpc form a complete set in 
space S. Another way to prove these results is to notice that the operator L is Hermitian 
and semi-negative definite by inspection and using the Cauchy-Schwartz inequality. Thus 
operator — 7U is Hermitian and semi-negative definite in the space S. Hence according 
to the Sylvester inertia law (see, e.g.. Theorems 4.5.8 and 7.6.3 in [20]), all eigenvalues of 
Cpc are real and non-positive, and the eigenfunctions of Cpc form a complete set in space 
S. As a result, under the time-step restriction At < Atmax, the convergence condition 
|1 + ApcAi| < 1 holds for all non-zero eigenvalues Apc of Cpc- On the other hand, if 
At > Atmax, 1 + ^min^t < — 1, hence iterations diverge. 

To complete the proof of Theorem 3, it remains to consider the kernel of jCpc in space S 
and verify that functions in this kernel do not affect the convergence of the iterations. 
For these functions, we have 

L^r - 7u = 0. (12) 
Taking the inner product between this equation and and noticing ^' € S", we get 

(L^',*)=0. (13) 

Since L is semi-negative definite, from the Cauchy-Schwartz inequality and the condition 
for its equality to hold, we get 

Li* = Pu. (14) 
On the other hand, differentiating Eq. (1) with respect to fi, we find that 

Li% = u. (15) 

Thus the solution \1' of Eq. (14) is equal to /?u^ plus functions in the kernel of Li. Due to 
the constraint e S and recalling our assumptions of (Uj^,u) = ^-P'(ju) ^ and u being 
orthogonal to the kernel of Li, we get (3 = 0. Thus the kernel of jCpc in space S is the 



20 



same as the kernel of Li. Since this kernel only contains invariance modes by Assumption 
1, it does not affect convergence of the iterations. Combining all the results obtained above, 
Theorem 3 is then proved. 



4.2 The power-conserving squeired-operator method for K equations with 
K propagation constants 



Next, we consider another class of equations (1) in the form 

Lqu = Loou - diag(/ii, . . . , jix) u = 0, 



(16) 



where u(x) = [1*1,^2, . . . , ukV is a real vector solitary wave, the superscript "T" represents 
the transpose of a vector, and are real propagation constants which are assumed to be 
independent of each other. The key feature of this case is that the number of independent 
propagation constants jik is equal to the number of components in the vector solution u(x) . 
Defining the powers of individual components (x) of the solution as 



and introducing the following K x K matrix 



D 



KxK 



1< k< 



(17) 



(18) 



then the PCSOM for Eq. (16) we propose is 

Pk 



{Uk 



Uk,n+1, ^<k<K, 



where 



u„+i = u„-M Uh\M %u-diag(7i,...,7K)u| At, 

>- J U=Un, IJ-k=IJ-k,n 

_ {uk, [M'^L|M-^Lou]fc) _ {uk, [M-iLoou]fc) 



(19) 



(20) 



(21) 



u=u„ 



Here [-jfc represents the fc-th component of the vector inside the bracket, \;\k,n represents 
the n-th iteration of the fc-th component of the vector inside the bracket, and M = 
diag(Mi, M2, . . . , Mk) is a positive definite and Hermitian acceleration operator. 
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The convergence properties of this PCSOM for Eq. (16) are similar to those of the previous 
PCSOM, and they are summarized in the following theorem. 



Theorem 4 Let all fi/. 's be independent, Assumption 1 he valid, u(x) he orthogonal to the 
kernel o/Li, det{D) ^ 0, and define ^tmax by Eq. (8), where A^j„ here is the minimum 
eigenvalue of the operator Cpcv defined as 

Cpcv^ = [tv^ - diagiii, 7]^)«} , (22) 

where 

Ly* = -L\m-^ [Li* - diag {Pi,..., Pk) u] , (23) 

and 

{uk, M^^Uk) ' {uk, M,:'^Uk) 

Then when At < Atmax> PCSOM (19)-(21) is guaranteed to converge to the solitary wave 
w(x) of Eq. (16) with component powers being [Pi, ... , PkY" > if the initial condition is close 
to w(x). When At > Atmax, PCSOM (19)-(21) diverges. 



Proof. The proof of this theorem is analogous to that of Theorem 3. By the linearization 
analysis, we find that the error of the iterated function satisfies the equation 

Un+l = Un + -CpcV Un Ai, (25) 

where operator Cpcv is as defined in Eq. (22). The power normalization (19) implies that 
it suffices to consider the eigenfunctions ^ of Cpcv satisfying the following orthogonality 
conditions: 

(^'fe,ufc)=0, k = l,...,K. (26) 

Similar to the operator Cpc in the proof of Theorem 3, we can show that in the space of 
functions satisfying the above orthogonality conditions, all eigenvalues of Cpcv are real and 
non-positive, and all its eigenfunctions form a complete set. The main difference between 
the previous operator Cpc and Cpcv here is in their kernels. The kernel of operator Cpcv 
contains functions * where Cpcv^ = 0- Similar to the case in Theorem 3, we can show 
that 

Li* = diag(/?i,/32,...,/?x)u. (27) 
Differentiating Eq. (16) with respect to ^i^, we get 

Liu^, = (0,...,ufe,...,0)^, k = l,...,K. (28) 
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Hence the solution of Eq. (27) is 

K 

* = ^/?,u^, (29) 

k=l 

plus the homogeneous solution of Eq. (27). Substituting this equation into the orthogonality 
conditions (26), and recalling the assumptions of Theorem 4, we get 

D • /3 = 0, (30) 

where matrix D is defined in Eq. (18), and (5 = {f3i,P2, ■ ■ ■ ■.PkY' ■ Since det(D) 7^ by 
assumption, the solution to Eq. (30) is /? = (J. Hence the kernel of Lpcv satisfying the 
orthogonality conditions (26) is the same as the kernel of Li. Thus, Theorem 4 is proved. 

The PCSOMs (3)-(5) and (19)-(21) have been presented above as separate cases because 
they have simple forms and also apply to many equations that arise frequently in ap- 
plications. For example, the method of Section 4.1 applies to linearly coupled nonlinear 
Schrodinger equations, while the method of Section 4.2 applies to nonlinearly coupled ones; 
sec also Examples 6.1-6.3 below. A PCSOM would also result for a more general case where 
in Eq. (1) of the form (16), the propagation constants /i^'s can be separated into several 
independent groups of equal /Xjfc's. An example of such a system of equations can be found, 
e.g., in [21]. In this case, the PCSOM can be constructed by combining the two PCSOMs 
(3)-(5) and (19)-(21). Here we group together the solution components n^'s whose corre- 
sponding propagation constants /i^'s are the same. These groups form sub- vectors in the 
whole solution vector u(x). Then the PCSOM for this more general equation is analogous 
to PCSOM (19)-(21), except that is replaced by each sub-vector of the solution, and 
Pfc replaced by that sub-vector's power. The convergence properties of this more general 
PCSOM are similar to those in Theorems 3 and 4, and will not be elaborated here. 

4.3 The squared-operator method with general quadratic conserved quan- 
tities 

In this subsection, we present the PCSOM for the case where the conserved quantities 
are not restricted to simple powers (as in Sees. 4.1 and 4.2). Rather, they can be more 
general quadratic quantities of the solutions (i.e. linear combinations of powers of individual 
components of the vector solitary wave). This case includes as particular cases the PCSOMs 
presented in Sections 4.1 and 4.2. We chose to present it separately from those two methods 
because its form is more complex than theirs, while the situations where neither of those 
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two simpler methods could apply are relatively rare. One example of such a situation is 
the system of coherent three- wave interactions in a quadratic medium [22, 23], which has 
only two conserved quantities which are linear combinations of the three wave powers (these 
conservation laws are known as Manley-Rowe relations). 



The class of problems described above can be formulated mathematically as follows: 

SQ 



Lqu = Loou 



Su 



11 = 0, 



(31) 



and 



Here u(x) 



K 



1=1 



j = l,...,r. (32) 

[ui, U2, . . . , is a vector solitary wave, 5/Su = [5/Sui, . . . ,5/5uk]'^ is the 



functional derivative, jl = (//i, . . . ,Hr)^ is the vector of r linearly independent propagation 
constants (1 < r < ii'), Pfe are defined in (17), qjk are constants, and Cj are the speci- 
fied values of the quadratic conserved quantities Qj. All variables and constants involved 
are real-vahicd. Without any loss of generality, we assume that the r X K matrix of the 
coefficients qj^ is in the reduced echelon form such that 



Qij 



0, i>j, 

1, i = j- 



(33) 



Then, using the terminology of textbooks in introductory Linear Algebra, we refer to powers 
{Pk}k=r+i {Pk}k=i as "independent" and "dependent" powers, respectively. 



The PCSOM that we propose for Eqs. (31)- (32) has the following form: 



Uk,n+1 



Ck-E 



K 

l=r+l 



QklPl 



Uk,n+1, 



{Uk,n+lT'^k,n+l) 

Un+1 = Un- M-^ fL|M-^Lou - B7„ 



l<k<r, 
At, 



(34) 



(35) 



U=U„, /l=/i„ 

where the notations for Li, M, and the subscripts are the same as before, Pk = {uk,n+i,Uk,n+i) , 



B 



5Q 

6u-' 



% = (B, M-1b)-1(B,M-1l1m-1Lou) 



(36) 
(37) 
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and 

/!„= (B, M-iB)-i(B,M-iLoou) . (38) 

u=u„ 

The power-normalization step (34) is written with the account of the reduced echelon form 
(33) of the matrix [qjk)- Note that this step affects only the "dependent" powers, while the 
"independent" ones do not need to be normalized. 

Convergence conditions of the PCSOM (34)-(38) are similar to those of the PCSOMs de- 
scribed above. Introducing the notation 

then convergence conditions of this PCSOM is summarized in the following theorem. 



Theorem 5 Let Assumption 1 he valid, all columns of the matrix B in (36) he orthogonal 
to the kernel of Jji, and dct(D) ^ 0. Also, let At^ax be given hy Eq. (8) where A^m 
now is the minimum eigenvalue of the operator CpcG defined in Eq. (40) below. Then for 
At < Atjnax! the PCSOM (34)- (38) is guaranteed to converge to the solitary wave m(x) of 
Eq. (31) which satisfies the constraints (32), provided that the initial condition is close to 
m(x). When At > Atmax, the PCSOM (34)-(38) diverges. 



The proof of this theorem follows the lines of Theorems 3 and 4, thus we will only sketch it 
below, emphasizing the differences from the proofs of the latter theorems. 

Proof. Linearizing Eq. (35), and noticing that the power normalization step (34) does not 
affect the linearized equations (as in Sections 4.1 and 4.2), we find that the error satisfies 
the iteration equation similar to (6), except that the iteration operator now becomes 

CPCG"^ = M'^ (Lpcg^ - B(B, M-iB)-i (B, M-^Lpcg*)) , (40) 

where 

tpcG"^ = -L|m-i (Li^' - B(B, M-1B)-1(B, M-^Li*)) . (41) 

The eigenfunctions ^' satisfy the orthogonality relation 

{Bj,^) = 0, j = l,...,r, (42) 

with Bj being the j-th column of B. These conditions follow from the quantities Qj being 
conserved and from the relation (36) between B and Q/s. 
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The operator LpcG is Hermitian. In addition, using the generahzed Cauchy-Schwartz 
inequahty for any matrix functions of Fi and F2: 

(Fi,Fi) > (Fi,F2)(F2,F2)-^(F2,Fi), (43) 

one can verify that Ij pcG is semi-negative definite. Then following the lines of the proof 
of Theorem 3, one then finds that in the space of functions satisfying the orthogonality 
conditions (42), all eigenvalues of CpcG are real and non-positive, and all its eigenfunctions 
form a complete set. 

Then it remains to consider the kernel of CpcG^ which satisfies the equation 

LpcG* - B(B, M-^B)-i (B, M-^LpcG*) = 0, (44) 

and show that the eigenfunctions of this kernel can only be those in the kernel of Li and 
thus would not affect the convergence of this PCSOM in view of our assumptions. To that 
end, we take the inner product between Eq. (44) and use the orthogonality relations 
(42), then recall that operator LpcG is semi-negative definite, and finally notice that the 
Cauchy-Schwartz inequality (43) becomes an equality if and only if Fi and F2 are linearly 
related by a constant matrix. This yields 

Li* = B^, (45) 

where /3 = (/3i, . . . , /J^)"^ is a constant vector. Noticing the relations obtained by differenti- 
ating Eq. (31) with respect to /Xj (1 < j < r), we see, similarly to (29), that the solution to 
Eq. (45) is 

* = E/?.%. (46) 

plus functions in the kernel of Li. Substituting this solution into (42) and recalling that 
by the assumptions of Theorem 5, columns of B are orthogonal to the kernel of Li and 
det(D) 7^ 0, we find that /? = 0. Thus the kernel of CpcG is the same as that of Li. 
Summarizing all these results. Theorem 5 is then proved. 

4.4 The squared-operator method for general equations with arbitrary 
conserved quantities 

The most general case is probably that the equations depend on the propagation constants 
/Xfe's in a general (but linear) way, and the specified conserved quantities of the solutions 
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are not restricted to powers or linear combinations of powers, but can be arbitrary. These 
equations and constraints can be written in the general form 

Lou = Loou-B(u)/x = 0, (47) 

and 

Qj{u) = Cj, j = l,...,r, (48) 

where /7, = {fii, . . . , fif)'^ is the vector of all independent propagation constants, B is a 
general matrix function of u, and Qj{u) are arbitrary functionals which are pre-specified. 
All quantities involved are real-valued. This system (47)-(48) generalizes the system (31)- 
(32) in the previous subsection in two significant ways. First, the functionals Qj are not 
restricted to the quadratic form (32) of u, but are allowed to be arbitrary functionals. 
For instance, one can seek a solution with a prescribed value of the Hamiltonian. Second, 
matrix B is not restricted to the special form (36) as in Eq. (31), but is allowed to be 
arbitrary functions of u. In this general case, counterparts of the power normalization step 
(34) become impossible. Thus, for the solution to have the pre-specified quantities (48), 
new ideas are necessary. Our idea is to replace the power normalization step with adding 
new terms into the iteration step (35) in such a way that when iterations converge, the final 
solution is guaranteed to meet the constraints (48). Our proposed scheme for this general 
case with arbitrary conserved quantities, which we denote QCSOM, is 

U„+i = Un - M 

where fin is defined by Eq. (38), and /i > is a user-specified free scheme parameter whose 
purpose will be explained after the proof of Theorem 6. The idea behind this scheme is that 
instead of minimizing the functional appearing on the right-hand side of Eq. (6) in Section 2, 
one minimizes a modified functional which equals the one from Eq. (6) plus additional terms 
5 Z^[Qj(u) — Cj]^. The acceleration (with the operator M"^) of this scheme is performed 
in the same way as the acceleration of (5). 

On the convergence of this QCSOM for Eqs. (47) and (48), we have the following theorem, 
which is very similar to Theorem 5 of the previous subsection. 

Theorem 6 Let Assumption 1 be valid, SQj/6u {j = 1, . . . r) be orthogonal to the kernel of 
Li, det(D) 7^ 0, and define At^ax by Eq. (8), where ^rnin here is the minimum eigenvalue 
of the operator Cqc defined in Eq. (50) below, then when At < Atmax, QCSOM (49) is 



lIm-^Lou + hY, (Qi(u) - Cj) 



6Qj 
Su 



At, (49) 



27 



guaranteed to converge to the solitary wave u(x) of Eq. (47) which satisfies the constraints 
(48), if the initial condition is close to m(x). When At > /^tmax> QCSOM (49) diverges. 

Here det(D) is as defined in Eq. (39), but Qj is an arbitrary function of u now. 

The proof of this theorem is also very similar to that of Theorem 5, thus we will only high- 
light its differences from the proof of the latter Theorem. Namely, the principal difference 
is that due to the absence of the normalization step in QCSOM (49), (38), there appears 
to be no counterpart of the orthogonality condition (42), which played a critical role in the 
proof in Section 4.3. However, as we will see, a direct counterpart of that condition will 
arise from different considerations. 

Proof. Linearizing the QCSOM, we find that the error satisfies the iteration equation 
similar to (6), except that the iteration operator now becomes 

Cqc^ = (^pcG^ - hp{^. , (50) 

where LpcG is defined by Eq. (41) of Section 4.3. The operator Cqc can be rewritten 
as Cqc = 'Mt^/'^Cqch M^/^. It is easy to check that Cqch is Hermitian. Using the 
generalized Cauchy-Schwartz inequality (43) , we can also verify that Cqch is semi-negative 
definite. Thus all eigenvalues of Cqc are real and non-positive, and all eigenfunctions of 
Cqc form a complete set. The kernel of Cqc satisfies the equation 

LpcG*-/.E(^,^)^ = 0. (51) 

Notice that both terms in the above equation are Hermitian and semi-negative definite 
operators, thus when taking the inner product between this equation and ^, we find that 
^ in the kernel of Cqc satisfies Eq. (45) as well as the orthogonality relations 

^)=0, i = l,...,r. (52) 

These relations are the counterparts of the orthogonality relations (42) of the previous 
subsection. Then the proof of this theorem is completed in exactly the same way as the 
proof of Theorem 5. 

Now we explain the reason for introducing the free parameter h into the scheme (49). Our 
experience shows that in some cases (especially when the conserved quantities Cj are large) , 
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the 5Qj/5\i terms in Eq. (50) with h = 1 cause operator Cqc's minimum eigenvalue A„j„ 
to be large negative — much larger than that of the operator — M^^Lgc" hi magnitude. 
This forces us to take very small time steps At (see Theorem 6), which severely slows down 
the convergence speed. When this happens, our strategy is to introduce a small parameter 
h into the scheme (49). The idea is to reduce the effect of the 6Qj/6u terms on the operator 
Cqc, and make its minimum eigenvalue close (or equal) to that of the operator — M^^Lqc". 
Hence the fast convergence of the scheme will be restored. This parameter h needs to be 
positive so that the relevant terms in Eq. (50) are semi-negative definite. 

To conclude this section, we point out that the general structures of the PCSOM (34)-(38) 
and QCSOM (49) can also be used to construct the imaginary-time evolution methods for 
the case when, in the above notations, 1 < r < K . To our knowledge, such a case was not 
considered in Refs. [9, 10, 12], where imaginary-time evolution methods for vector equations 
were reported. Of course, the corresponding imaginary-time evolution methods, unlike the 
PCSOM (34)-(38) and QCSOM (49), will not be universally-convergent. 



5 Squared-operator iteration methods for isolated solitary 
waves 

In many dissipative wave systems such as the Ginsburg-Landau equation, solitary waves 
exist only when the propagation constants in the equation take discrete (isolated) values. 
We call the corresponding solutions isolated solitary waves. For these isolated solutions, the 
propagation constants or powers of the solutions are unknown and need to be computed 
together with the solitary wave, thus the numerical schemes discussed in previous sections 
do not apply. In this section, we propose squared-operator iteration methods for isolated 
solitary waves. 

For simplicity, we consider the case where a vector isolated solitary wave exists when a 
single propagation constant takes on discrete values (other cases can be easily extended). 
The equation for the solitary wave can be written as (1), but the solution u(x) now exists 
only at an unknown discrete fj, value. We propose the following squared-operator iteration 
method for these isolated solitary waves (SOMI): 



Un+l = U„ — 



M-^l|m-^Lou1 At, (1) 



= /^n + (u,M-iLou) At, (2) 

U=U„, IJ.=IJ-n 



29 



where Li is the hnear operator as defined in Eq. (2), and M is a positive-definite and 
Hermitian acceleration operator. Linearizing this method around the isolated solitary wave, 
we get the iteration equation for the error (u„, as 



where 



]=il + AtCj) , (3) 

( Un\( -M-^LIM-I (LiU„ - flnU) \ , . 

' [fin )-[ (u, M-i(Lifi„-/l„u)) )■ ^ ^ 

It is easy to check that operator Cj can be written as diag(M"i/^ 1) Cjh diag(M^/^ 1), 
where Cjh is Hermitian and semi-negative definite. Thus all eigenvalues of Cj are real 
and non-positive, and all its eigenfunctions form a complete set. The kernel of Cj contains 
functions [<I>(x),0]"^ with $(x) in the kernel of Li, as well as functions [F(x),l]^ where 
F(x) satisfies the equation 

LiF = u, (5) 

and must be bounded. Assuming that u is not orthogonal to the kernel of operator L| 
(which is the generic case), then Eq. (5) has no bounded solution, hence the kernel of £/ 
only contains the invariance modes [$(x),0]"^. Then under Assumption 1, SOMI (l)-(2) 
will converge if At < —2/Ajnin and diverge if At > —2/Amin, where A^j^ is the minimum 
eigenvalue of operator Cj. 

Following the motivation for MSOM (1) in Sec. 3, we can construct the modified squared 
operator method for isolated solitons in Eq. (1) as follows (MSOMI): 



and 
where 



Un+l — Un + 



— Mn + 



-M-^lIm-ILou - anOnGn] At, (6) 

{u,M-^LoVL) - anenHn] At, (7) 

J U=U„, /i=Mn 



1 1 

" {MGn,Gn) + m ~ ((LiG„ - HnU) , M'^ (LiG„ - HnU))At' ^ ^ 

On = -(LiG„ - Hnu, M-^Lou), (9) 

and {Gn,Hn) are functions specified by the user. A good choice, which is a counterpart of 
Eq. (5), is 

Gn = e„ = U„ — U„_i, Hji = — fJ-n-l- (10) 

We can show that under similar conditions as in Theorem 2 and the assumption below Eq. 
(5), this MSOMI also converges for all isolated solitary waves. 
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6 Examples of applications of squared-operator iteration meth- 
ods 



In this section, we consider various examples of physical interest to illustrate the perfor- 
mances of the proposed schemes. Most of these examples are two-dimensional, thus the 

shooting method can not work. For some of the examples such as Examples 1 (b,c) and 5, 
other numerical methods such as the Petviashvili method and the imaginary-time evolution 
method can not work either. 

Example 6.1 Consider solitary waves in the two-dimensional NLS equation with a periodic 
potential, 

U:,^ + Uyy - Vo (sin^ X + sin^ y) U + a\U\'^U = -fiU, (1) 

which arises in nonlinear light propagation in photonic lattices and Bose-Einstein condensa- 
tion in optical lattices. Here /it is the real propagation constant, and U is & complex solution 

in general. Writing U into its real and imaginary components, U{x,y) = u(x,y) + iv{x,y), 
we get equations for the real functions u and v as 

Uxx + %j/ ~ Vo (sin^ X + sin^ y^ u + a{u^ + v^)u = —jiu, (2) 

Vxx + Vyy — Vq {sit? X + sin^ y^ v + a{u'^ + v'^)v = —jiv. (3) 
This system is rotationally invariant, i.e., if {u,v)'^ is a solution, so is 

cos 9 — sin 6 \ I u \ 
sin 9 cos 9 I \ ^ I ^ 

where 9 is the angle of rotation (which is constant). This invariance, in Eq. (1), corresponds 
toU —>■ Ue^^. This rotational invariance induces an eigenmode {—v,u)'^ in the kernel of Li. 
Clearly, this eigenmode is orthogonal to the solution (n, u)^, satisfying the orthogonality 
condition in Theorem 3. The kernel of Li docs not contain other eigenfunctions (at least 
in the generic case) since there are no other invariances in this system. 

Solitary waves in Eq. (1) exist only inside the bandgaps of the system. Thus, the bandgap 
information is needed first. The bandgap diagram at various values of Vq is displayed in 
Fig. 2(a). For illustration purpose, we fix Vq = 6, and determine the solitary waves in 
different bandgaps below. 
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(a) Vortex solitons in the semi-infinite bandgap under focusing noniinearity: 

For focusing noniinearity, £7 = 1. In this case, Eq. (1) admits various types of real and 
complex solitary- wave solutions in every bandgap of the system [24, 25, 26]. Here we 
determine a vortex-soliton solution at // = 3 (P = 14.6004), which is marked by letter 
'a' in the semi-infinite bandgap of Fig. 2(a). This solution is complex valued with 
a non-trivial phase structure, and it is displayed in Fig. 2(c, d). Similar solutions 
have been reported theoretically in [24] before, and have since been experimentally 
observed [27, 28]. To determine this solution, we apply the SOM (4), MSOM (1), (5), 
and PCSOM (3)-(5) on Eqs. (2)-(3), starting from the initial condition 

(4) 

In addition, we choose the acceleration operator M as 

M = (c-5,,,.-9^,)diag(l,l). (5) 

The spatial derivatives as well as are computed by the discrete Fourier transform 
(i.e., by the pseudo-spectral method). The computational domain is — 67r < x,y < Gtt, 
discretized in each dimension by 256 grid points. It should be noted that the size of 
the computational domain and the number of grid points have very little effect on 
the convergence speed; they mainly affect the spatial accuracy of the solution. For 
these three schemes, we found that the optimal (or nearly optimal) parameters are 
(c. At) = (3.7, 0.8) for SOM and PCSOM, and (c. At) = (3.8, 0.6) for MSOM. At 
these choices of parameters, the error diagrams versus the number of iterations are 
displayed in Fig. 2(b). Here the error is dofinod as the difference between successive 
iteration functions: 

en = ^J {Un - Un-1, Un - Un-l) ■ (6) 

We see that all three schemes converge rather quickly. The convergence speeds of 
SOM and PCSOM are almost the same, but MSOM converges much faster. It should 
be noted that the amount of computations in one iteration is different for these three 
methods, with the ratio roughly of 1 : 1.7 : 2 for SOM, PCSOM, and MSOM. When 
this factor is also considered, we conclude MSOM converges the fastest, with SOM 
second, and PCSOM third. 

(b) Solitons in the first bandgap under defocusing noniinearity: 

Next, we consider solutions in the first bandgap (between the first and second Bloch 
bands) under defocusing noniinearity {a = —1). For this purpose, we pick /x = 5, 
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marked by letter 'b' in the first bandgap of Fig. 2(a). At this point, Eq. (1) admits 
a real- valued gap soliton, which is displayed in Fig. 3(a). Similar solutions have been 
reported in [29, 30] before. To determine this solution, we apply the SOM (4), MSOM 
(1), (5), and PCSOM (3)-(5) on Eq. (1), starting from the initial condition 

U{x, y) = 1.15 sech(x^ + y^) cos(x) cos(y). (7) 

We take M as 

M = c - dxx - dyy. (8) 

The computational domain is — Stt < x, y < Stt, discretized in each dimension by 128 
grid points. For these three schemes, we found that the optimal (or nearly optimal) 
parameters are (c,At) = (1.8,0.6) for SOM and PCSOM, and (c,At) = (2.9,1.7) 
for MSOM. At these choices of parameters, the error diagrams versus the number of 
iterations are displayed in Fig. 3(b). Similar to the vortex soliton in Fig. 2, we find 
that the convergence speeds of SOM and PCSOM are almost the same, but MSOM 
converges much faster. 

(c) Vortex solitons in the second bandgap under defocusing nonlinear ity: 

We now determine vortex solitons in the second bandgap under defocusing nonlinearity 

(cr = —1). For this purpose, we pick fi = 9.4, marked by letter 'c' in the second 
bandgap of Fig. 2(a). At this point, a vortex soliton with distinctive amplitude and 
phase distributions exists (see Fig. 4(a, b)). It is noted that this vortex soliton is not 
of the type reported in [30, 26], which lie in the first bandgap (in our notations). To 
our knowledge, this type of soliton has never been reported before in the literature. To 
determine this solution, we apply the SOM (4), MSOM (1), (5), and PCSOM (3)-(5) 
on Eqs. (2)- (3). The initial condition is the vortex solution of these equations at a 
different value = 9.42, which we in turn obtained by the continuation method from 
a small-amplitude solution near the edge of the bandgap. The acceleration operator 
M is the same as (5). The computational domain is — IOtt < x,y < lOvr, discretized 
in each dimension by 256 grid points. For these three schemes, we found that the 
optimal (or nearly optimal) parameters are (c. At) = (4.2, 1.7) for SOM and PCSOM, 
and {c,At) = (4,3.1) for MSOM. At these choices of parameters, the error diagrams 
versus the number of iterations are displayed in Fig. 4(b). In this case, MSOM is 
again the fastest. However, unlike the above two cases, PCSOM converges much 
slower than SOM now. 

Example 6.2 Consider the following system arising from a second-harmonic generation 
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1 2 

VXX + ^Vyy +2'" 



liv. (10) 



(SHG) model, 

UXX + Uyy +UV = liU, (9) 

1 

2 

Solutions of similar systems have been considered in a number of studies; see, e.g., a recent 
paper [31] and references therein. Note that in the original reduction from the SHG model, 
the right hand side of the u-equation is usually 2/xu. But in order to cast the equations 
into the form (1), we have divided the v-equation by 2, so that we work with Eq. (10) 
instead. Here we deliberately make the Vxx and Vyy coefficients different, so that the radial 
symmetry is broken, hence this system is not reducible to a radially symmetric (and hence 
essentially one-dimensional) problem. At /Lt = 0.1, this system admits a solitary wave with 
total power P = {u,u) + {v,v) = 47.3744, which is displayed in Fig. 5(a, b). We take the 
initial condition as 



u{x,y) = sech.{0. 35 \J x'^ + y"^), f = 0.3sech(0.5y + y^). (11) 
The acceleration operator M is taken as 



M = diag 



A* dxx dyy , /X + ^ ^ dxx 2 ^yy 



(12) 



where the choice of constants /x and /x + ^ is motivated by our earlier results on optimal 
accelerations for imaginary-time evolution methods [13]. The computational domain is 
—25 < x,y,< 25, and the number of grid points along each dimension is 64. The iteration 
results of SOM (4), MSOM (1), (5) and PCSOM (3)-(5) at optimal (or nearly optimal) At 
values 0.37, 0.59, and 0.63 respectively are displayed in Fig. 5(c). Again, MSOM delivers 
the best performance, and PCSOM is the slowest. 

Example 6.3 The next example is the coupled two-dimensional NLS equations with sat- 
urable nonlinearity, 

Uxx + Uyy + — r^" — oT^ = /^i"' (1^) 

which has been studied, e.g., in [9, 32, 33]. Here s is the saturation constant. At s = 
0.5, /Lt = 1,^2 = 0.5 (Pi = 85.3884,^2 = 29.1751), this system admits a solution whose 
n-component is single-humped, but its v-component is a dipole state. This solution is 
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displayed in Fig. 6(a, b). Note that this solution is not radially symmetric, thus is not 
reducible to a one-dimensional problem. Taking the initial condition as 

u{x, y) = 3e-°-2^' , v{x, y) = 1.5r e-''-^''" cos 9, (15) 

where (r, 0) are the polar coordinates, the acceleration operator M as 

M = diag [ill - dxx - dyy, /X2 - dxx - dyy] , (16) 

the computational domain as — 12 < x, y, < 12, and the number of grid points along the two 
dimensions as 64, the iteration results of SOM (4), MSOM (1), (5) and PCSOM (19)-(21) 
at optimal At values 1.9, 2.65, and 1.85 are displayed in Fig. 6(c). As in (12) above, the 
choice of the constants and /i2 in the acceleration operator (16) is also motivated by our 
previous studies on the accelerated imaginary-time evolution method [13]. Again, MSOM 
delivers the best performance, and PCSOM is the slowest. 

Example 6.4 The next example is intended to compare the performances of PCSOM (34)- 
(38) and QCSOM (49). This example comes from the three-wave interaction system [22, 23] 
and has the form 

Uxx + Uyy + VW = fJ,iU, (17) 
Vxx + Vyy + UW = fJ,2V, (18) 
Wxx + Wyy +UV = {m + H2)w, (19) 

where u, v and w are real functions, and /xi and iJ,2 are propagation constants. At /xi = 0.5 
and /X2 = 1, this system has a radially symmetric solution displayed in Fig. 7(a), and 

Qi = {u, u) + {w, w) = 66.3096, Q2 = {v, v) + {w, w) = 47.2667. (20) 

Here we want to determine this solution with the pre-specified quantities Qi and Q2 as 
above. Note that this problem is of the form (31)-(32) and (47)-(48), but not of the form 
(1) or (16). Taking the initial condition as 

it(a;, y) = 2.5sech0.8r, i;(,x, y) = 2.2sech0.8r, w{x, y) = 1.9sech0.8r, (21) 

where r = vx^ + y^, the acceleration operator M as 

M = diag [m - dxx - dyy, H2 - dxx - dyy, 111 + 112- dxx - dyy] , (22) 

the computational domain as —15 < x,y,< 15, the number of grid points along the two 
dimensions as 64, and the parameter h in the QCSOM as h = 0.01, the iteration results of 
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PCSOM (34)-(38) and QCSOM (49) at the (same) optimal value At = 0.49 are displayed 
in Fig. 7(b). This figure shows that the QCSOM converges slightly slower than the PC- 
SOM. However, it is noted that each QCSOM iteration involves less computations than the 
PCSOM, thus we conclude that the performances of PCSOM and QCSOM are comparable. 

Example 6.5 The Ginzburg-Landau equation is of the form 

i^t + (1 - lii)^xx - ho^ + = 0, (23) 

where 70 is the damping/pumping coefficient (when 70 is negative/positive). We seek 
solitary waves in this equation of the form 

^x,t) = U{x)e''^\ (24) 

where /x is a real propagation constant, then function U{x) satisfies the equation 

{l-jii)U^^-i^oU + \U\^U = fiU. (25) 

If 70 or 71 is non-zero, then solitary waves in this equation are always complex- valued, and 
they can only exist at isolated propagation constants. When 70 = 0.3 and 71 = 1, the 
solitary wave, which exists at the discrete value /U = 1.2369, is plotted in Fig. 8(a). Writing 
this equation into two real-valued equations for [Re(C/), Im(C/)], and applying SOMI (l)-(2) 
or MSOMI (6)-(10) with M = c — dxx and initial conditions uo{x) = 1.6sech(x), no = 1.2, 
we can obtain this isolated solitary wave. At optimal scheme parameters (c. At) = (1.6, 0.3) 
for SOMI and (c, At) = (1.4, 0.12) for MSOMI, the error diagram is displayed in Fig. 8(b). 
Here the error is defined as 

en = \J {Un - Un-l-Ahi " «ri-l) + \^^n - l^n-l\- 

We see that MSOMI converges much faster than SOMI, which is consistent with previous 
numerical experiments. 

Wc should point out that in this example, since 70 > 0, the soliton we obtained is unstable 
(because the zero background is unstable). So this soliton can not be obtained by simulating 
the Ginzburg-Landau equation (23). However, our proposed method SOMI/MSOMI can 
produce this unstable solution quite easily. 



7 Summary 

In this paper, we have developed three iteration methods — the squared operator method 
(SOM), the modified squared operator method (MSOM), and the power (or any quantity)- 
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conserving squared operator method (PCSOM/QCSOM), for computing solitary waves in 
general nonlinear wave equations. The solitary waves can exist at either continuous or 
discrete propagation constants. These numerical methods are based on iterating new differ- 
ential equations whose linearization operators are squares of those for the original equations. 
We proved that all these methods are guaranteed to converge to all types of solitary waves 
as long as the time step in the iteration schemes is below a certain threshold value. Due to 
the use of acceleration techniques, these methods are fast converging. Since these methods 
are compatible with the pseudo-spectral method, their spatial accuracy is exponentially 
high. Furthermore, these methods can treat problems in arbitrary dimensions with little 
change in the programming, and they are very easy to implement. To test the relative 
performances of these methods, we have applied them to various solitary wave problems 
of physical interest, such as higher- gap vortex solitons in the two-dimensional nonlinear 
Schrodinger equations with periodic potentials and isolated solitons in Ginzburg-Landau 
equations. We found that MSOM delivers the best performance among all the methods 
proposed. 

Even though MSOM delivers the best performance, SOM and PCSOM/QCSOM have their 
own advantages as well. For instance, SOM is simpler to implement. PCSOM/QCSOM 
would be advantageous if the problem at hand specifies the power or other conserved quan- 
tity of the solution rather than the propagation constants. In addition, PCSOM/QCSOM 
works for linear eigenvalue problems (by setting the linear eigenfunction to have a fixed 
norm), while SOM and MSOM do not. In some cases, SOM or PCSOM is more tolerant to 
the choice of initial conditions, i.e., they converge for a larger range of initial conditions than 
MSOM. Thus all the methods developed in this paper can be useful for different purposes, 
and the reader can pick them judiciously depending on the problem at hand. 

It is noted that these methods can deliver good performance even with suboptimal choices 
of scheme parameters. But how to get the best performance out of these schemes (i.e., how 
to find optimal scheme parameters) is still an important open question. For the Petviashvili 
method, conditions for optimal or nearly optimal convergence have been studied in [6, 14]. 
For the accelerated imaginary-time evolution method, optimal acceleration parameters have 
been obtained for a large class of equations [13]. Such results can help us to select the scheme 
parameters for the squared-operator methods in this paper (as we have done in Examples 
6.2-6.4). But a rigorous and comprehensive study on optimal scheme parameters for the 
squared-operator methods proposed in this paper is a non-trivial issue and will be left for 
future studies. 
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Appendix: Families of squared-operator methods 

The squared operator methods proposed in this paper were based on the operator L given 
in Eq. (9). It turns out that one can construct a family of squared operator methods which 
contain the methods in this paper as particular cases. Consider the following squared 
operator iteration method for Eq. (1): 



Un-t-l = Uri 



M-"LtM-^Lou] At, (A.l) 
^ Ju=u„ ^ ' 



where M is a positive definite Hermitian operator, and a and h are arbitrary constants. 
The linearized equation of this method for the error u„ is still Eq. (11), except that L is 
replaced by 

Lf = -M-"l|m-''Li (A.2) 
now. This operator can be rewritten as 

Cf = -M-"/^ (M-^'/^LiM-^/^y (M-''/2LiM-"/2j M"/^, (A.3) 

thus its eigenvalues arc clearly all non-positive. Repeating the proof in Sec. 2, we can 
readily show that the SOM (A.l) is guaranteed to converge if At < —2/Amin, where A^m 
is the minimum eigenvalue of operator Cf. If we choose a = 1 and 6 = 1, then method (A.l) 
becomes SOM (4). For the family of SOMs (A.l), the corresponding MSOMs, PCSOMs, as 
well as methods for isolated solitary waves can be readily constructed. Below we illustrate 
such a construction for a particular choice of a and b. Namely, we consider the member of 
family (A.l) with a = and b = 2, whose implementation of each iteration requires slightly 
fewer operations than the implementation of method (4) in the main text. In this case, the 
corresponding squared operator methods are listed below. 

• SOM for Eq. (1): 

Un+i = Un - [lIm-^LouI At. (A.4) 



u=u„ 



38 



MSOM for Eq. (1): 

Un+l = U„ 

where 



lIm-^Lou - an{Gn, l1m-2Lou)G„ 

1 1 



u=u„ 



At, 



a„. = 



(G„,G„) (M-iLiG„, M-iLiGn)Ai' 
and Gn is a user-specified function such as (5). 

PCSOM for Eq. (1): 



(A.5) 
(A.6) 



Un+l = U„ — 



(Un+1, U„+i) 



u 



n+l) 



and 



u=u„, ;i=/i„ 

2- 



At, 



(u, l1m-2Lou) (u, M-"Loou 

7= -/ ^ ) ^ 



(u, u) 



(u, M-2u) 



(A.7) 
(A.8) 

(A.9) 



We note that for this power-conserving scheme, the 7 term in Eq. (A.8) can be 
dropped due to the presence of the power normahzation step (A.7). It is easy to check 
that the reduced scheme has the same hnearized iteration operator as the original one 
above, thus possesses the same convergence properties. However, this can not be done 
for the PCSOM (3)- (5) in the main text if M / 1. The other forms of the PCSOM 
for the cases (16) and (31)-(32) can be similarly written down. 



QCSOM for Eqs. (47)- (48): 



Un+l — U^ 

where 



lIm-^Lou + hj2 {Qj{u) - C,) ^ 



At, (A.IO) 



U=U„, IJ.=IJ.„ 



/xn= (B, M-2b)-1(B,M-2Loou) 
SOMI for isolated solitary waves in Eq. (1): 



Un+l = u„ 



l|m-2LoU 



U=U„, /U=A(„ 



At, 



Un+l = lJin+ (u, M ^Lqu) 



U=U„, /U=/i„ 



At, 



(A.ll) 

(A.12) 
(A.13) 
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MSOMI for isolated solitary waves in Eq. (1): 

u„+i = u„ + f-L|M-2LoU - a„0„GJ At, (A.14) 



— /^n + 

where 



(u,M-^Lou)-a„6'„i/„ At, (A.15) 

JU=U„,;i=/i„ 



_ 1 1 

" (G„, G„) + H^~ (M-i (LiG„ - Hnu) , M'^ (LiG„ - HnU))At' ^ 

On = -{Gn, l1m-2Lou) + (F„u, M-^Lqu) , (A.17) 
and {Gn,Hn) are user specified functions such as (10). 

All these methods can be shown to have similar convergence properties as their counterparts 
in the main text. 

As we have already noted, in all the methods presented in this Appendix starting with 
Eq. (A. 4), each iteration involves a little less computations than their counterparts in the 
methods presented in the main text. However, we did not advocate for these methods for 
two reasons. First, we can show that the convergence of these methods is always slower 
than that of their counterparts in the text when Li is Hermitian. In such a case, using 
Theorem 5.6.9 in [20], we can show that ^min/^max is smaller for C than it is for £j, 
hence SOM (4) converges faster than its counterpart (A. 4) in view of Eq. (15). Second, 
for the squared-operator methods in the main text, we can carry out explicit convergence 
analysis on some familiar examples such as the NLS equation. This would not be possible 
for the methods considered in this Appendix. Overall, our numerical testing shows that the 
squared-operator methods presented in the main text give the best performance among the 
family of methods (A.l) with other choices of a and b. 
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Figure 1: Analysis of the SOM (4) applied to the soliton (17) of the NLS equation (16). 
(a) Eigenvalues of operator M~^Li, with M given in (18), versus the acceleration pa- 
rameter c; (b) maximum and minimum (nonzero) eigenvalues of the iteration operator 
C = — {M~^Li) ; (c) graph of convergence factor R^{c) versus c; its minimum occurs at 
c = Copt = 6 — vTS; (d) convergence factor function R{At; c) versus At at c = Capt- 
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Figure 2: (a) Bandgap structure of Eq. (1); (b) error diagrams for SOM (4), MSOM (1), 
(5) and PCSOM (3)-(5) at optimal c and At values (see text); SOM and PCSOM are 
almost indistinguishable; (c, d) a vortex soliton in the semi-infinite bandgap of Eq. (1) with 
focusing nonlinearity. Here /i = 3 (P = 14.6004), marked by letter 'a' in panel (a), (c) is 
the amplitude {\U\) distribution, and (d) the phase distribution. 
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Figure 3: (a) A solitary wave in the first bandgap of Eq. (1) with defocusing nonlinearity 
at /i = 5 (P = 2.4936), marked by letter 'b' in Fig. 2(a); (b) error diagrams for SOM 
(4), MSOM (1), (5) and PCSOM (3)-(5) at optimal c and At values (see text); SOM and 
PCSOM are almost indistinguishable. 
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Figure 4: (a, b) A vortex-soliton solution in the second bandgap of Eq. (1) with defocusing 
nonhnearity at ^ = 9.4 {P = 18.3578), marked by letter 'c' in Fig. 2(a). (a): amplitude 
distribution; (b) phase distribution, (c) error diagrams for SOM (4), MSOM (1), (5) and 
PCSOM (3)-(5) at optimal c and At values (see text). 
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Figure 5: (a, b) A fundamental solitary wave in the SHG model (9)-(10) at fi = 0.1 
[P = 47.3744); (c) error diagrams for SOM (4), MSOM (1), (5) and PCSOM (3)-(5) at 
optimal scheme parameters (see text). 
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Figure 6: (a, b) A dipole-mode vector solitary wave in the coupled NLS system (13)-(14) at 
m = 1 and = 0.5 (Pi = 85.3884, P2 = 29.1751); (c) error diagrams for SOM (4), MSOM 
(1), (5) and PCSOM (19)-(21) at optimal scheme parameters (see text). 
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Figure 7: (a, b) A fundamental soliton in the three-wave model (17)-(19) at jii = 0.5 and 
112 = 1, where Qi = {u,u) + {w,w) = 66.3096, and Q2 = {v,v) + {w,w) = 47.2667; (c) 
error diagrams for PCSOM (34)-(38) and QCSOM (49) at optimal scheme parameters (see 
text) . 
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Figure 8: (a) An isolated solitary wave in the Ginzburg-Landau equation (25) with 70 = 0.3 
and 71 = 1; (b) error diagrams of SOMI (l)-(2) and MSOMI (6)-(10) at optimal c and At 
values (see text). 
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